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Abstract 

Although stream-dwelling gudgeons (Cyprinidae, genus: Gobio) are widespread in Central Europe, the taxonomy of this 
group and the distribution of its species are still unexplored in detail. The aims of our study are to ascertain taxonomic 
composition and distribution of the former Gobio gobio superspecies in the inner area of the Carpathian Basin. Since the 
presence of cryptic species is suspected in this area, we examined the taxonomic and phylogenetic relationships of Central 
European Gobio taxa by sequencing the mitochondrial DNA control region (mtCR). Additionally, we characterized the 
genetic structure of 27 stream-dwelling gudgeon populations of this area by Amplified Fragment Length Polymorphism 
(AFLP). Results of mtCR analysis proved the presence of three species already known as G. obtusirostris (dominant in NW- 
Hungary), G. gobio (sporadic) and G. carpatliicus (sporadic). Additionally, the analysis revealed the existence of one doubtful 
taxon, G. sp7 (dominant in NE-Hungary), and a new isolated haplogroup (dominant in SW-Hungary). Although Network 
analysis showed significant detachment among haplogroups, their genetic distances were quite small. Therefore Bayesian 
phylogenetic analysis showed weak nodal support for the branching pattern both for newly described haplotypes, and for 
the already accepted species. AFLP data showed distinct population structure and a clear pattern of isolation was revealed 
by distance of stocks. At the same time, level of separation was not affected by the altitudinal position of sites. Moreover we 
found three major clusters of populations which were separated according to hydrographic regions, and corresponded to 
the findings of mtCR analysis. Our results suggest the on-going speciation of gudgeons in the Carpathian Basin, however 
the separation of haplogroups seems to only be an intermediate phase. The discovered natural pattern seems to be only 
slightly influenced by anthropogenic impacts. Additionally our results put into question the suitability of the recently 
accepted within Gobio genus taxonomy. 
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Introduction 

Freshwaters are exceedingly diverse ecosystems, but at the same 
time they are extremely sensitive to habitat degradation and 
pollution [1]. Accurately quantifying their taxonomic and 
functional diversity is a fundamental requirement of conservation 
biological, ecological, biogeographical, and macroevolutionary 
research [2], [3], [4]. However, when species are not clearly 
distinguishable by the conventional methods using ecological and 
morphological traits, have highly similar environmental needs and 
reveal a high level of phenotypic plasticity, then estimates of 
biodiversity and ecosystem functioning will be biased [5] , [6] , [7] , 
[8]. 

With the increasing use of molecular techniques, it has become 
evident that many species formerly believed to have widespread 
geographical distribution can in fact be divided into numerous 
more or less discrete entities -so called cryptic or sibling species by 
Mayr [9]- or represent genetic gradients between separating 



species (Le. on-going speciation) [10], [11], [12], [13], [14]. Such 
phenomena are more likely to occur in organisms with limited 
dispersion ability and/or in organisms living in separated or 
narrowly connected habitats [15], [16], [17]. Many stream- 
dwelling fishes have specific environmental needs and therefore 
form discrete populations, not only between geographical areas 
with separated catchment systems, but also between closely related 
sites of the same catchment [18], [19], [20], [21]. These isolated 
fish populations may then genetically differentiate with time, 
although they may still maintain their similar morphological 
appearance and ecological function [22]. However, this process is 
not yet fully understood in seemingly well connected catchment 
systems. 

The type species of the Gobioninae subfamily (Fam: Cyprinidae), 
the common gudgeon Gobio gobio Linnaeus (1758), was known as 
the most widely distributed lentic gudgeon species in West Eurasia 
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[23]. However, the high between and within population morpho- 
logical variability of this taxon [24] resulted a long-standing debate 
regarding its taxonomical status [25], [26]. As a result, numerous 
forms/varietas/.subspecies were described, nearly all from larger 
catchment areas within its range (e.g. [23]). Recent genetic studies 
have raised some of the former sul)spc(:ies to species level and 
several new species have been described as well [27], [28] [29], 
[30], [31], [32]. However, the taxonomic and genetic status of 
gudgeon populations in the Carpathian Basin is still unclear. Gohio 
gohio was considered a common species in the waters of the 
Carj)atluan Basin for a long time [33], [34]. On the contrary, 
recent studies [31], [35] excluded the Carpathian Basin from the 
potential range of G. gobio and suggested the occurrence of the 
Danube gudgeon, Gobio obtusirostris Valenciennes (1842) in the 
western region of the basin, and the Carpathian gudgeon, Gobio 
carpathicus Vladykov (1925) in the drainage system of River Tisza. 
Furthermore, Mendel et al. [31] indicated the presence of a 
"species" (sic!: Gobio spl) from the catchment of River Tisza, but 
the taxonomic position of this newly described taxon has not yet 
been firmly established. It is also important to note that the 
findings of Kottelat and Freyhof [3,5] are based exclusively on data 
from the literature. Moreover, although the study of Mendel et al. 
[31] is the most comprehensive genetic study on the Middle 
European Gobio species to date, it included only a very limited 
number of samples from the Carpathian Basin and furthermore all 
of those samples originated from the edges of this region. The 
above mentioned Gobio taxa (i.e. G. gobio, G. obtusirostris, G. 
carpathicus and G. spl) are morphologically very similar [35], [36] 
and thus their distribution and ecology- cannot be explored without 
molecular justification [24]. Furthermore, recent studies on G. 
gobio and the related species [37], [38], [39], [40] found 
remarkably high genetic and morphologic variability between 
and within catchment areas, supporting the likelihood of the 
presence of cryptic species. 

This study aims to ascertain taxonomic composition and 
distribution of the former G. gobio superspecies in the inner 
Carpathian Basin. Specifically, we (i) examine the taxonomic and 
phylogenetic relationships of Gohio taxa and the presence of cryptic 
species, and (ii) characterize the genetic structure of gudgeon 
populations with special attention paid to the effects of hydrolog- 
ical distance and elevation as possible forces facilitating genetic 
separation. 

To unravel the taxonomic relationships of stream-dwelling 
gudgeons inhabiting the central area of the Carpathian Basin we 

use the same methodology and molecular marker (sequencing the 
mitochondrial Control region) as was used by Mendel et al. [31], 
thus our results are comparable with their findings. Moreover, we 
also screened for Amplified Fragment Length Polymorphisms 
(AFLP's) as a supplementary method of analysis for genome-wide 
genetic variation [41]. 

Materials and Methods 

Ethics Statement 

This study was carried out following relevant national and 
international guidelines pertaining to the care and welfare of fish. 
Collections were made by electrofishing, partially from sampling 
sites which are located within protected areas. Moreover, each 
species within the Gobioninae subfamily is protected in Hungary. 
Electrofishing in protected areas and any procedure (collection 
and storage of tissue samples) to be applied to protected species are 
subject to authorisation in Hungary. Fin tissue collection and 
storage were approved by the National Inspectorate for Environ- 
ment, Nature and Water, Hungary (permission numbers: 14/ 



3714-2/2009, 14/1237/2/2010, 14/881/5/2011, 14/678-9/ 
2012). Fish collected for this study were narcotized using clove 
oil. After fin tissue samphng, when they regained consciousness, 
they were returned to the wild. Field studies did not involve fish 
that were endangered (The lUCN Red List of Threatened Species 
V. 2013.1; www.iucnredlist.org). 

Study Area 

All sampling sites were situated in the inner area of the 
Carpathian Basin, which belongs to the drainage system of the 
Danube River. Based on its topographic characteristics, the 
Hungarian part of Middle Danubian hydrosystem can be divided 
into two larger catchments and ten smaller sub-catchments 
(Table 1, Figures 1 and 2). The hydrography of this area is 
biaxial. From the western region, all watercourses empty into the 
Danulx; Ri\'cr. The (eastern part of the country belongs to the 
drainage system of the Tisza River, which is the largest tributary of 
the Danube (157 000 km^ catchment area). The structure of this 
drainage system is dendritic, with the Tisza River forming the 
central axis. All the studied watercourses connected to the middle 
section of the River Tisza [42], therefore this region was not 
differentiated further. The hydrography of the Danubian system is 
more complicated. This system consists of three comparatively 
isolated subsystems: North, Middle and South Danubian regions 
(Table 1). North Danubian region is formed by the drainages of 
River Raba (Raab), River Ipoly (Ipel) and by the drainages of 
some direct inflowing streams. Middle Danubian catchment 
originally joined to the River Danube through a marshy area. 
Until the construction of the Sio canal at the end of the 19th 
century, there was only intermittent connection between the 
Danube and the Lake Balaton drainage system. Therefore this 
subdrainage had been hydrologically isolated to some degree from 
the others until the last century. Waters from South Danubian 
region flow into the River Drava (Drau), which empties to the 
Danube River at Osijek (Croatia). 

Fish Sampling 

Fish samples were colletled between 2009 and 2012 by 
elertrofishing from 27 sites across five sub-i:atchments of the 
Danube River and five sub-tributaries of the Tisza River (Table 1, 
Fig. 1). Since individuals of the Gobio genus show notable 
phenotypic plasticity, we investigated only adult (>60 mm 
standard length) specimens characterised by Gobio gobio- like 
morphological features, such as dispersedly spotted dorsal and 
caudal fins, and with no epithelial crests on the scales situated on 
the predorsal region of the body [35]. 

Molecular Methods 

DNA extraction and purificatian. Fin clips of 241 speci- 
mens were sampled and stored in 96% ethanol at — 20°C until 
DNA extraction. DNA was isolated with a DNeasy Blood and 
Tissue kit (Qiagen, Germany), using 10-20 mg of fin tissue as per 
the manufacturer's instructions. QuEility and quantity of the 
extracted DNA were verified using a NanoDrop 2000c Spectro- 
photometer (Thermo Scientific, USA). 

Mitochondrial sequence data. DNA of 168 out of 241 
individuals (111 from the Danube River and 57 from the Tisza 
River catchments) were used for the amplification of the 
mitochondrial control region (mtCR). The sequences of mtCR 
were amplified by polymerase chain reaction (PCR) using the 
primers CR159 (CCCAAAGCAAGTACTAACGTC) and 
CR851 (TGCGATGGCTAACTCATAC) ([3,3]). PCRs were 
carried out using 0.2 \x\ of 5 U/|J,1 Taq DNA polymerase 
(Fermentas), 2.5 ^ll of lOX Taq buffer, 1.7 ^ll MgClj (25 mM), 
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Figure 1. Location of tlie Carpathian Basin in Europe (A), location of the sampling sites (1-27) in the Carpathian Basin (B) and PCoA 
representation of hydrologic distances between the sampling sites (C). Hungarian country border is marked with the dotted line. Different 
symbols refer to sites belonging to different catchment areas: O - North Danubian; □ - Middle Danubian; O - South Danubian and A - the Tisza River 
catchment. For detailed information see Table 1. 
doi:1 0.1 371 /journal.pone.0097278.g001 



0.2 |il dNTPs (10 niM), 0.3 [il of each primer (20 [iMj, 2.0 |il 
template DNA, and 17.8 |J,1 purified and distilled water in a final 
volume of 25 |xl. Reactions were performed in a MJ Research 
PTC-200 Peltier Thermal Cycler under the following conditions: 
95°C for 1 min, followed by 37 cycles of 94°C for 45 s, annealing 
at 52°C for 30 s, and an extension temperature of 72°C for 45 s, 
followed by a final extension at 72°C for 8 min. PGR products 
were purified from 1 % agarose gel using the MUlipore Ultrafree- 
DA DNA extraction kit. PGR products were sequenced on an ABI 
3730XL sequencing machine by MWG-Biotech AG (http://www. 
mwg-biotech.com). Sequences were edited manually and aligned 
using the program Geneious 5.4 [43] and ClustaLX 2.0.11 [44]. 
Newly described haplotypes have been deposited in the Gene- 
Bank. Calculation of sequence polymorphism and haplotype 
detachment was performed using DnaSP 5.10 software [45]. 
Sequence divergence was calculated with net nucleotide diver- 
gence (Da) in MEGA5 [46]. 

AFLP. To verify the results of the mitochondrial CR 
sequencing, a complementary method, Amplified Fragment 
Length Polymorphism (AFLP), was carried out; which is a 
reproducible, PCR-based molecular genetic method [41]. Alto- 
gether, 241 specimens were surveyed according to the following 
protocol. 200 ng DNA extracted from caudal fm tissue was 
digested at 37°C for 2 hours in a final volume of 10 |iL with 2.5 U 
Msel, 5 U EcoRI enzymes and 2 |J.L NEBuffer4 (New England 
BioLabs, USA). Enzymes were then inactivated at 65°G for 
20 min. Adaptor ligation was carried out at 24°G for 16 hours in 
20 |J,L final volume containing the total digestion mixture, 



0.25 nM EcoRI, 2.5 |xM Msel adaptors, 200 cohesive end unit 
T4 Ligase and 1 x T4 Ligase Reaction Buffer (New England 
BioLabs). After heat inactivation at 65°C for 10 min, 10 |J.L of 
digested, hgated mixture was diluted 10 fold with nuclease free 
water. Pre-selective PGR was carried out with AmpliTaqGold 360 
Master Mix (Applied Biosystems, USA) in a final volume of 20 |J.L 
containing 0.5 \iM Eco-A (5' GAGTGCGTACCAATTGA 3'), 
0.5 nM Mse-G primer (5' GATGAGTCCTGAGTAAG 3') and 
5 (iL diluted, digested, ligated DNA. The PGR was started at 
94°G for 2 min followed by 20 cycles of 94°G for 30 sec, 56°G for 
1 min, 72°G for 1 min and a final elongation at 72°G for 7 min. 
Selective PGR was performed in a final volume of 20 |J.L 
containing AmpliTaqGold 360 Master Mix, 0. 1 |xM fluorescendy 
labeUed Eco-AGT primer (5' 6FAM GAGTGGGTAGGAATT- 
GAGT 3'), 0.25 ^M Mse-GTT primer (5' GATGAGTGGT- 
GAGTAACTT 3'), 2 ixL PGR product from the pre-selective 
PGR. Gycling conditions of the touchdown PGR were as follows; 
enzyme activation at 94°G for 2 min followed by 13 cycles for 
30 sec at 94°G, for 30 sec at 65°G and decreased by 0.7°G in each 
cycle, and for 1 min at 72°C, then 23 cycles for 30 sec at 94°G, for 
30 sec at 56°G and for 1 min at 72°G, followed by 5 min at 72°G. 
Digestion, ligation and PGRs were carried out in a Gene Amp 
PGR System 9700 (Applied Biosystems). Fragment analysis was 
performed on an ABI 3 1 30 sequencer (Apphed Biosystems, USA) 
and data were analysed with Peakscanner vl.O (Applied Biosys- 
tems). Electropherograms were automatically analysed with 
tinyFLP [47], by the following scoring parameters: min. height: 
90, max. width: 1, min. size: 50, max. size: 500, range (+/~): 0.5, 
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Figure 2. Median-Joining network of mtCR sequence data relating Gobio spp. with previously published data. Circle size is relative to 
the number of individuals carrying the same haplotype. Line length refers to the genetic distances of haplotypes. Small open circles represent median 
vectors (missing or theoretical haplotypes). CR01-1 7: Haplotypes of the 168 specimens analysed in this study. Letter code of haplogroups/"boxes" 
and numbers (N°) of previously published haplotypes in each box correspond with the numbers and codes displayed in Table 2. 
doi:1 0.1 371 /journal.pone.0097278.g002 



min peak-peak dist.: 1, peak height diflFerence: 0, min. freq.: 0.1, 
max. freq.: 90. From the 553 peaks detected in total, 154 selected 
bands were retained. After further evaluation (e.g. specimens with 
a small number of peaks were excluded from the analysis) a dataset 
of 229 specimens was used for further statistical analyses. 

Data Analysis 

Mitochondrial sequence data. To shed light on the 
taxonomic relationships, alignment of all haplotypes found in this 
study and the previously pubhshed Gobio haplotypes (source: [31]) 
described from the neighbouring regions (e.g. Central Europe, 
Balkan Peninsula, and Anatolia) was performed (Table 2). 
Originally the sequences revealed in this study were 65 1 bp long, 
but for the Network analyses we had to align them to the GenBank 
sequences of these closely relative Gobio species. Therefore for the 
Network analyses a 652 bp long dataset were used. Network was 
constructed using the median-joining algorithm in Network v. 4.6. 
[48]. Similar haplotypes were classified arbitrarily into hap- 
logroups (see "boxes" in Fig. 2). Dilferentiation within and among 
haplogroups was tested by analysis of molecular variance 
(AMOVA; [49]) with 9999 permutations. 

To construct the phylogenetic tree, the (652 bp long) sequence 
set analyzed in the Network analysis was aligned against further 
haplotypes used as outgroup taxa of varying putative phylogenetic 
depths (sources: [50], [51], [31], and Mendel et al. unpublished 
data): Romanogobio vladykovi (GenBank a.n.: EF427385), Romanogobio 
banaticus (GenBank a.n.: EF427393), Sarcocheilichthys variegatus 
(GenBank a.n.: NC004694), Rhodeus ocellatus kurumeus (GenBank 
a.n.: AB070205). Thus the lenght of aligned sequence set was 
666 bp. Bayesian phylogenetic analysis was conducted by Markov 



chain Monte Carlo method (B/MCMC), and it was performed in 
MrBayes 3.2 [52]. The best-fitting models of DNA substitution 
were selected for analysis using Akaike information criterion (AIC) 
implemented in the jModelTest 0.1.1 [53], [54]. jModelTest 
indicated that Hasegawa, Kishino and Yano substitution model 
[55] with gamma-distributed rate heterogeneity (a = 0.57 10) 
(HKY-hG) was the best fitting. We conducted Bayesian tree 
construction with 6 chains, 2 independent runs and 7 million 
generations. Trees were sampled every 1000th generation. The 
first 1 0000 generations were discarded as burn-in. We plotted the 
log-likelihood scores of sample points against generation time 
using Tracer 1.5 [56] to ensure that stationariness was achieved 
after the first 10000 generations by checking whether the log- 
likelihood values of the sample points reached a stable equilibrium 
plateau. We used the remaining trees with average branch lengths 
to create a 50% majority-rule consensus tree with the sumt option 
of MrBayes. Posterior probabilities were obtained for each clade. 

AFLP analysis. Higher level differentiation of Gobio assem- 
blages was assessed using STRUCTURE 2.3.3 [57] to estimate 
the most probable number of genetic groups (clusters, K) for all 
analysed individuals. Values of K were investigated from 1 to 10, 
with a burn-in period of 10000 followed by 100,000 iterations and 
10 runs for each K using an admixture model with correlated 
allele frequencies. Results of the Bayesian statistics were evaluated 
by Structure Harvester [58] implementing the Evanno method 
[59]. To characterise the standard measures of population genetic 
diversity, the percentage of polymorphic loci (%), mean unbiased 
heterozygosity, and unbiased Nei's gene diversity [60] were 
calculated. Within population genetic distance (GD) was calculat- 
ed using the following equation: 
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Figure 3. Bayesian consensus tree derived from the analysis of the mtCR sequence data. Haplotypes revealed in this study are marked 
with their CR codes (see Table 3). Bayesian posterior probabilities are assigned on nodes. For GenBank accession numbers see the text and Table 2. #: 
taxon name described in Mendel et al. (2008). 
doi:1 0.1 371 /journal.pone.0097278.g003 



Table 2. GenBank haplotypes used for the network computatii 
displayed in Figure 2. 



where 2n^y equals the number of shared character states and n is 
the total number of binary characters. Population genetic structure 
was characterised by hierarchical AMOVA [49] with 9,999 
randomisations. Isolation by distance was estimated by a Mantel 

1. Numbers (N°) of haplotypes correspond with the numbers 





GenBank accession numbers 


code 


taxon name by GenBanit 






N-3 


N-4 


N-5 


'A' 


Gobio obtusirostris 


EU131554 


EU131557 


EU131558 






'C 


Gobio spl^ 


EU131564 


EU131565 


EU131563 






'D' 


Gobio gobio 


EU131542 


EU131544 


EU131543 


EU131545 


EU131546 


'E' 


Gobio skadarensis 


EU131568 


EU131569 


EU131567 






'F' 


Gobio carpathicus 


EU131561 


EU131552 


EU131560 


EU131559 




'G' 


Gobio ohridar)us 


EU131572 


EU131571 


EU131573 


EU131570 




'H' 


Gobio insuyanus 


EU131576 


EU131574 


EU131578 


EU131580 


EU131579 



(#: taxon name described in Mendel et a!., 2008). 
doi:10.1371/journal.pone.0097278.t002 
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Figure 4. Determination of tKie number of clusters best fitting the AFLP data: STRUCTURE-based mean±SD liiceiihood values of ten 
runs for each K from K=1 to 10 (A), similarity coefficient (min., mean and max values) of ten runs for each K from K=1 to 10 
(mean±SD) (B), Delta K statistic (C), and Triangle plot with allocation of individuals to clusters mapped according to K = 3 (D). Where 
O: specimens from the North Danubian region (the two emphasized individuals are identified as Gobio gobio in the CR analysis), □: Middle Danubian 
sites, O: South Danubian sites, A: specimens from the Tisza region (the two emphasized individuals are identified as Gobio carpathicus in the CR 
analysis). 

doi:1 0.1 371 /journal.pone.0097278.g004 
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Figure 5. Plots of cross-validation tables for AFLP data based 
on DAPC. Correctly classified individuals are placed on the diagonal. 
The square size equals the number of individuals of posterior group 
assignment based on posterior probabilities. Rows correspond to actual 
sites (a priori), while columns correspond to inferred sites (posteriori). 
Squares with broken lines show regional detachments of populations. 
doi:10.1371/journal.pone.0097278.g005 



test [61] using pair-wise Opt (similar to pairwise Fj,, [62] data). Nei 
unbiased genetic distances and pairwise hydrological distances 
were derived from a hydrological map (1:10 000) with 9,999 
randomisations. All of these calculations were made in GenAlEx 
v6.5 [63] statistical software. The inbreeding coefficient (FJ and 
frxation index, as a measure of poptdation dilferentiation (Fs,) from 
AFLP markers, were computed using the Bayesian ABC4F 
software [64]. The percentage of polymorphic loci, mean unbiased 
heterzygosity, F;s and F^, values of the studied populations (where 
N&8) were compared with the altitude of the sampling sites by 
Spearman rank correlations. We calculated the membership 
probabilities of each individual for the different a priori groups 
(i.e. populations, where N^8 and geographical regions) based on 
retained discriminant functions using cross-validation with Dis- 
criminant Analysis of Principal Components (DAPC) [65]. AH the 
files used for statistical analyses are available in the supplementary 
material. 

Results 

Mitochondrial DNA Sequence Analysis 

Aligned sequences of 651 3'- end CR mtDNA were obtained 
from 1 68 individuals grouped into 1 7 haplotypes. Sequence data 
of 1 5 previously undescribed haplotypes are highhghted in bold in 
Table 3. These have been deposited in the GenBank database 
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Figure 6. Recent river network of the Danubian region (Hungary). The four subcatchment areas are indicated by different colours. Flow 
direction of Zaia River during the Middle and Late Pleistocene is indicated by blackframed green and red arrows respectively. The thick broken line 
indicates the Pleistocene watershed between the North and IVliddle Danubian regions. Intra-valley drainage divides (Sikhegyi, 2002) are shown by 
black bidirectional arrows. Numbered circles: sampling sites displayed in Table 1. Dotted arrows: recent flow direction. 
doi:1 0.1 371 /journal.pone.0097278.g006 



under Accession Nos. KC757328-42. The sequences of the CR03 
and CR14 haplotypes had already been identified from the 
Danube catchment in the Czech Republic and in Slovakia [33], 
and were demonstrated by megaBLAST [66] to be Gobio gohio 
(100% similarity with the specimen: EU131542) and G. carpathicm 
(100% similarity with the specimen: EU131559) respectively. 
From the North Danubian region only three and from the Tisza 
River catchment a total of five haplotypes were displayed. The 



Middle and South Danubian regions were the richest in 
haplotypes. Moreover nine out of the 10 haplotypes found were 
unique to these regions (Table 3). Results of the median-joining 
Network analysis showed that haplotypes described from Hungary 
were classified into five haplogroups (A, B, C, D, F 'boxes' in 
Fig. 2). Through the AMOVA analysis, 85% of the total genetic 
variance was explained by among haplogroup differences, and in 
addition significant (p<0.001) differentiation was found in each 
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pairwise comparison of haplogroups, confirming the arbitrary 
classification pattern. 

With the exception of haplogroup 'B', each one corresponds to 
a previously described "species" (Table 4). In the Network 
analysis, 46 investigated specimens were identified as G. ohtusirostris 
(haplogroup 'A'), 55 as G. spl described by Mendel et al. [31] 
(haplogroup 'C'), two as G. gohio (haplogroup 'D') and two as G. 
carpathicus (haplogroup 'F') (Fig. 2, Table 4). Altogether 63 
specimens originating from the Middle and South Danubian 
region form a distinct, currently unidentified group (haplogroup 
'B'), which was in a transitional position between G. spl and G. 
ohtusirostris (Fig. 2). The results of Bayesian phylogenetic analysis 
(Fig. 3) were similar to those obtained through the Network 
analysis. However, the posterior probabilities of nodes showed a 
high level of uncertainly (weakly supported branching) at both 
"lower" (i.e. species) and at "higher" (i.e. genus) levels as well. 

AFLP. The mean of the estimated Ln probability values from 
STRUCTURE analysis of the final matrix strongly increased 
between K = 1 and K = 3 and then consolidated at higher K values 
(Fig. 4a). The comparative statistics [67], [68] supported three 
major clusters (Fig. 4b, c). A triangle plot of the results (Fig. 4d) 
showed that individuals classified into the Cluster 1 originated only 
from the Danubian region while the individuals in Cluster 2 
originated completely from the Tiszanian catchment area. The 
Cluster 3, consisting of mainly Middle and South Danubian 
specimens, shows a transitional position (continuous transition) 
between the two aforementioned clusters. This pattern was similar 
to the hydrological distances between the sampling sites (see 
Fig. Ic). 

For population genetic analyses, data were used of those 21 
populations (196 individuals) where N&8. The average number of 
bands per specimen was 40.5±9.4 (ranging between 16.0 and 
62.0). Base population genetic features as: P. loci %, UHe, Fj^, F^t 
are given for each population where N^B in Table 3. Kruskall- 
Wallis tests revealed that the North, Middle and South Danubian 
and Tiszanian groups of populations did not show any significant 
differentiation in terms of their population genetic features 
(Table 3). 

Within population genetic distance (CD), pairwise C>pt and 
pairwise Unbiased Nei Genetic Distances data are displayed in 
Table 5. GD ranged between 25.4 and 38.8 (av. ±SD = 34.8±3.1) 
and neither differed significandy by subregion, nor correlated with 
the altitudinal position of the collection site. Mean Unbiased Nei 
Genetic Distances ranged between 0.007 and 0.100 (av. 
±SD = 0.035±0.020). AMOVA analysis showed that among 
group differences accounted for 12°/i> of the total genetic variance, 
and 193 out of the 210 pairwise comparisons (93%) showed 
significant (p<0.05) differentiation. Pairwise Opt data ranged 
between zero and 0.321 (av. ±SD = 0.1 16±0.06). This markedly 
strong population separation was verified by the results of D AFC. 
Assessing the AFLP dataset, 86% and 99% of the individuals were 
grouped correctly on population and region levels respectively in 
multidimensional space based on the cross-validation procedure 
within DAPC (Fig. 5). Results of Spearman rank correlations 
supported that the population genetic variables (P. loci %, UHe, 
Fjs, and F^^ were not significandy correlated (p<0.05) with the 
altitudinal position of the sites. 

According to the results of Mantel tests, hydrological distances 
correlated significantly and positively both with Unbiased Nei 
genetic distances (Rxy = 0.499, p<0.01) and with the pairwise Opt 
values (Rxy = 0.551, p<0.01) of the selected (Ns8) populations. 



Discussion 

Our results only partially support the earlier assumptions of the 
taxonomic composition and distribution of Gobio species in the 
inner Carpathian Basin [3 1] , [35] , and show that Gobio reveal high 
diversity, both taxonomically and from the population genetic 
standpoint, within this relatively small and well-connected 
catchment area. At the same time, the presence of more cryptic 
species is indicated and the fine-scale separation of the identified 
genetic lineages between sub-catchments supports the existence of 
ecological barriers and on-going speciation in Hungarian drainage 
systems. 

Taxonomic and Phylogenetic Features 

All of the Gobio haplotypes found in the inner Carpathian basin 
can be classified into the north European clade and all except one 
haplotype (G. gobio) belong to the north-eastern European subclade 
described by Mendel et al. [31]. At the same time, our results 
indicate that the taxonomic status of the stream-dwelling gudgeons 
inhabiting the inner area of the Carpathian Basin is more complex 
than was previously presumed. Although no remarkable differ- 
ences were detected in the morphological and meristic traits of the 
specimens analyzed, haplotypes of three previously described 
species (G. obtusirostris, G. gobio, and G. carpathicus), a doubtful taxon 
(G. spl) and an additional, distinct haplogroup were distinguished 
from the study area. 

Although the haplotype of G. carpathicus occurred in the Tisza 
River catchment, this area was dominated by the haplotypes of G. 
spl (Fig. 2, Table 4). Gudgeon stocks inhabiting the Danubian part 
of the country showed greater taxonomic complexity. Contrary to 
the earlier hypotheses [35], we found the haplotype of G. gobio in 
the Carpathian Basin (in Cuhai-Bakony-er). Furthermore, haplo- 
types of G. obtusirostris proved to be dominant only in the North 
Danubian region. A distinct, but highly diverse haplogroup 
(haplogroup 'B') was dominant in the Middle and South Danubian 
regions (Fig. 2, Table 4). 

DifiFerentiation of haplogroups 'A' and 'B' (Table 4) may be 
attributed to a population split caused by paleohydrographic 
changes that took place in the geologic recent past. Namely, 
approximately 140,000 years ago a new watershed with an east- 
west direction formed, separating the North and Middle Danubian 
region [69] . This changed the flow direction of the Zala River, 
originally flowing northward following the current channel of the 
Marcal River, southward to the Drava River. Separation of the 
Middle and South Danubian regions began only at the end of the 
Pleistocene, 14-16,000 years ago by the formation of Lake Balaton 
[70], [71]. Namely, the Zala River, and all the smaller streams 
flowing southward until then, turned toward this newly formed 
depression. The phylogenetic effect of the relatively old watershed 
separating the North and Middle Danubian region and the 
incomplete splits - "intra-valley drainage divides"- (Fig. 6) between 
the Middle and the Southern regions [72] was proven by the fa[:t 
that (1) the Network analysis showed significant diflerentiation 
between haplogroup 'A' and 'B', but (2) haplotypes did not differ 
between the drainages of Lake Balaton, River Kapos-Sarviz- 
system, and River Drava (Fig. 2, Table 4). 

Although some specimens characterised by CROl (G. obtusiros- 
tris) haplotypes also occurred in the Middle Danubian region, their 
presence was restricted to only those sites that were in the vicinity 
of fish ponds and/or where trout {Oncorhynchus and/or Salmo) 
stocking occurred [73]. Therefore, we assume unintentional G. 
obtusirostris introductions with gamefish (trout) stocking and thus a 
secondary, anthropogenic contact between the phylogenetically 
separated haplogroups. 
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< 



Haplogroup 'B' showed similar genetic distances from the G. 
obtusiwstris (1 .3 1 % ± 0. 3 1 %) and from G. .f/^i (1 .6 1 % ± 0. 1 7 %), as an 
interspecific genetic divergence among some already accepted 
gudgeon species such as G. skadarensis and G. carpathicus 
(1.29% ±0.18%) or G. skadarensis and G. ohridams 
(1.38%±0.17%) (Table 4). These dilTerences make species level 
detachment of the haplogroup ^B' or G. spl taxon as well. At the 
same time, the results of Bayesian phylogenetic tree analysis 
showed weakly supported differentiation among haplogroup 'B' 
and some already accepted Gobio species (e.g. G. obtusirostris), and 
did not clearly support the recently accepted within-genus 
taxonomy. Our results indicate that the reproductive isolation of 
these entities may have only began in the geohistorical recent past, 
presumably in the middle Pleistocene. Contrary to the findings of 
the Network analysis, Bayesian phylogenetic computations in most 
cases query the species level differentiation within the Gobio genus. 
This premise is supported by the fact that in the case of other fish 
species, e.g. topmouth gudgeon - Pseudorasbora parva (Temminck & 
Schlegel, 1842), stone loach - Barbatula barbatula (Linnaeus, 1758), 
> and grayling - Thymallus thymallus (Linnaeus, 1758), a similar or 

S higher degree of differentiation among haplogroups is considered 

S to be not more than subspecies level detachment [20], [74], [75]. 

.1 In addition, some authors [76], [77] have suggested that the 

^ genetic distance of haplogroups must be greater than or equal to 

S ten times the level of within-haplogroup differences to distinguish 

"I separate species. In our study, the G. insuyanus is the only taxon 

which fulfils the above mentioned criteria (Table 4). 

I Population Genetic Variables 

g The values of basic population genetic parameters (P. loci%, 

!fl UHe, Fst, and GD) did not show significant differences among 

subdrainage basins. Moreover, none of these variables showed 
.§ significant correlation with the altitudinal position (as a possible 

g marker for population isolation due to differences in die habitat 

g use of fish) of the sampling sites. In the case of basic population 

genetic parameters, the local environmental conditions and the 
degree of hydrographic isolation are likely to be more important 
than either the altitudinal position or the taxonomic arrangement 
of the inhabiting specimens. 

The population genetic features of gudgeon assemblages 
inhabiting the northwest region of Hungary differed notably from 
the other studied Hungarian assemblages. It is the only area where 
I 5 statistical analyses suggested considerable gene flow (Fig. 5, 

y Table 5). This may be attributable to the species level differences. 

At the same time, landmarks of the river systems characterising 
I this area assure convenient migration routes among sites, 

c Similarly, the occurrence of Gobio gobio (haplotype) may be 

S attributed to the role of Danube River. At the other Danubian 

^ Regions the population structure was much more differentiated. 

For the Middle Danubian Region, notable differences were found 
among closely situated sites. Results of the mtCR analyses revealed 
I the existence of a different haplotype (CROl) in addition to the 

S g assumed "native" haplotype group from this area (Table 3). 

Therefore these differences may be caused by accidental G. 
obtusirostris introductions to this area. 

Results of STRUCTURE analysis are in accordance with the 
results of mtCR sequence analyses. Both inferred the existence of 
three larger clusters/haplogroups within the Carpathian Basin. 
Furthermore, both analyses indicated the transitional position of 
Cluster 2/haplogroup 'B' between the Cluster 1 /haplogroup 'A' 
and the Cluster 3/haplogroup 'C The two specimens identified 
as G. gobio and two specimens identified as G. carpathicus by mtCR 
analyses did not show notable separation from the others by AFLP 
analysis (Fig. 4D) which may suggest interspecific hybridisation in 
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these cases. There are numerous reports of interspecific and 
intergeneric [Romanogobio and Gobid) hybridisation [23], [78], [79] 
of European gudgeon species and our results support these 
findings. 

Mantel test results revealed a clear pattern of isolation by 
hydrographical distance. Taxonomic and population genetic 
differences of the studied Gobio stocks were simultaneously changed 
by the growing hydrographical distances. This natural pattern is 
just sUghtly diminished by anthropogenic impacts. 

Consequendy, the genetic analyses confirmed the results of 
former analyses, which were based on mainly morphologic/ 
morphometric variables [23], since they revealed that the Middle 
European Gobio "species" form an extremely diverse and variable 
group. At the same time explanation of the phylogenetic 
relationships and within-genus taxonomic features are still pardy 
unresolved. Our results showed that because of the casual 
immigration and/or the accidental introductions, and the 
sympatric occurrences, the location of the collection site is not a 
convenient feature to discriminate these "species" occurring in 
Hungarian waters. 

Our results indicate that these cryptic Gobio entities form a 
relatively young phylogenetic group and that the genetic 
differences among them are not strong enough in most cases for 
species level differentiation. Moreover, considering the possibilities 
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of interspecific and intergeneric hybridisation, the recent taxo- 
nomic partitioning of the Gohio genus needs re-evaluation. 
However, the discovered genetic diversity is probably very 
vulnerable. Since the separation of haplogroups seems to be only 
an intermediate phase of an on-going speciation and stream- 
dwelling gudgeons have specific environmental needs and a 
restricted habitat area at present, habitat alteration and accidental 
stocking may easily damage the integrity of haplogroups. 
Consequendy, conservation actions should be implemented to 
preserve the exceptional diversity of this fish group. 
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